function [Ytk1, Ytk2] = aggregate_sparseY(Y, T, K1, K2)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% This function takes a sparse vector Y of dimension T x (K1*K2) and takes the sum of Y along dimension k1 and k2.
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Inputs:
	% Y:			T x (K1*K2)  (sparse)
	% T:            integer
	% K1:            integer
	% K2:            integer
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Outputs:
	% Ytk1:         T x K1
	% Ytk2:         T x K2
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	%%% Compute Ytk1
	Y = reshape(Y,[T*K1 K2]);
	Ytk1 = reshape(sum(Y,2), [T K1]); % T x K1
	
	%%% Compute Ytk2
	Y = reshape(Y,[T*K1*K2 1]);
	idxes = find(Y);
	[tt_vec,k1_vec,k2_vec] = ind2sub([T K1 K2], idxes);
	tk2_vec = sub2ind([T K2], tt_vec, k2_vec);
	Ytk2 = accumarray(tk2_vec, Y(idxes), [T*K2 1]);
	Ytk2 = reshape(Ytk2,[T K2]);
end
